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Abstract. This paper studies the bond valence method (BVM) and its 
application in the non-isovalent semiconductor alloy (GaN)i_x(ZnO)x- Particular 
attention is paid to the role of short-range order (SRO). A physical interpretation 
based on atomic orbital interaction is proposed and examined by density- 
functional theory (DFT) calculations. Combining BVM with Monte-Carlo 
simulations and a DFT-based cluster expansion model, bond-length distributions 
and bond-angle variations are predicted. The correlation between bond valence 
and bond stiffness is also revealed. Finally the concept of bond valence is extended 
into the modelling of an atomistic potential. 
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1. Introduction 

Isovalent semiconductor alloys are extensively studied 
since their structural, electronic and optical proper¬ 
ties can be tuned by varying the alloy composition. 
More recently, non-isovalent semiconductor alloys be¬ 
gin to attract attention. For example, the pseudobi¬ 
nary (GaN)i_x(ZnO)x alloy is attractive for its high 
efficiency in photocatalytic water splitting[l]. From 
the theoretical perspective, my co-workers and I have 
recently predicted strong short-range order (SRO) in 
the (GaN)i_x(ZnO)x alloy due to its non-isovalent 
nature[2]. The role of SRO in determining the atomic, 
electronic and vibrational properties is also revealed. 
To fulfill the local charge neutrality, the substitutional 
SRO is accompanied with and compensated by atomic 
deviation from the ideal lattice positions. Therefore it 
is imperative to study the composition-, temperature- 
and SRO-dependent (a;,T, 11) structural relaxations. 
In Ref. [2] the (GaN)i_x(ZnO)x alloy is efficiently 
represented by a SRO-modified version of the Special 
Quasirandom Structure (SQS)[3] approach. The SQS- 
based approach “mimics” the statistics of the corre¬ 
lations. The present study aims directly at statisti¬ 
cally reliable predictions of bond-length distribution 
and bond-angle variation. 

The bond valence method (BVM) is widely adopted 
in solid state chemistry for various applications includ¬ 
ing prediction of molecular geometry [4], construction 
of atomic potentials for perovskite oxides [5] [6], and cal¬ 
culation of the acidity constant piFa[7][8]. Its power 
in predicting the energetics for non-isovalent semicon¬ 
ductor alloys is recently demonstrated[9] [10]. In inor¬ 
ganic chemistry the BVM is commonly recognized as 
an empirical tool, the underlying physics of which is 
not widely discussed. For example, the fact that bond 
valence correlates strongly with bond length[4] reflects 
the connection between bond valence and bond-length- 
dependent transferable force constant [11]. Also the 
correlation between total energy and bond valence [4] is 
not yet fully understood. Brown[4] proposed a “more 
rigorous but less physical” analogy of the Kirchhoff 
circuit law which treated the bond valence network 
as a capacitive electric circuit. Burdett[12] interpre- 
tated BVM from the molecular orbital basis. There is 
also some similarity between the bond valence and the 
Mayer bond order[13]. In the present study, particu¬ 
lar attention is paid to the theoretical standing of the 


BVM. Density-Functional Theory (DFT) calculations 
are performed wherever necessary. 

The present study also has a computational motiva¬ 
tion. Electronic structure methods can handle fairly 
large supercells (e.g., over 100 atoms). For non- 
isovalent semiconductor alloys, a large supercell is fa¬ 
vored in order to average out the fluctuation error due 
to the finite size. As large structural relaxations are 
expected, accurate DFT total energy and force calcu¬ 
lations are computationally expensive. Therefore it is 
desirable to pre-relax the internal atomic positions in 
an economical way. The BVM-predicted bond lengths 
and bond angles suit this purpose well. 


2. The Bond Valence Method 

The BVM is extensively discussed in Ref. [14]. 
Each nearest-neighbor cation-anion bond is assigned 
a bond valence vjj. Next nearest-neighbor cation- 
cation/anion-anion interactions are neglected. The 
bond valence sum (BVS) of an atom is defined as the 
sum of the bond valences surrounding the atom. Each 
atom has an ionic valence V equal to its corresponding 
formal ionic charge. By convention, V(Ga)=-|-3, 
V(N)=-3, V(Zn)=+2, V(0)=-2, vij = -vjj. Of 
crucial importance for non-isovalent semiconductor 
alloys are two rules: (1) the valence sum rule V{T) = 
'^jVij, and (2) the valence loop rule J2ioop''^iJ ~ 
0. The valence sum rule is an equivalent statement 
of the principle of local charge neutrality, with the 
correlation Pjj ex vij where Pjj is the Mulliken 
overlap population[15][16]. The valence loop rule is 
also known as the equal valence rule, since the zero 
circulation condition is equivalent to the minimization 
oiEpjvjj ( see for example the appendix in Ref. [10] ). 
The solution is a set of {vjj} which minimizes the 
measure of the total energy E = a j is 

the correlation constant) under the constraint of the 
valence sum rule. 


As for the measure of the total energy, in solid state 
language, a perturbation expansion of the orbital 
interaction energy reads 


iei iei 





( 1 ) 
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where e° and (/)° denote atomic energy and orbital 
respectively. Capital I,J and lowercase i,j refer to 
atomic and orbital indices respectively. Assuming 
the correlation Hij cx Sij {Hij is the matrix element 
{4>^\H 'and Sij is the overlap integral the 

relaxation energy E — '^ti is then approximately equal 
to ajj j^j ^iji where the denominator e° — e° is 
absorbed into a correlation constant ajj. The overlap 
integral (f)) can be expressed in a separable 

form Sij{r)f{9,(j))[12]. The angular dependence is 
lifted after summing J2iei jeJ 

interacting orbital pairs. The summation over orbital 
pairs then reduces to the summation over atom pairs. 
Finally the measure of the total energy E = a j vjj 
is obtained, with the assumptions vjj ~ Sjj = 

j^j ^ij — O!, while a is to be fitted 

by DFT total energy calculations. In the analogy of 
the Kirchhoff circuit law, the bond capacitances are all 
equal [4], which is equivalent to assuming a equal for 
different types of atomic pairs. 

The radial dependence of Sij{r) leads naturally to 
the empirical exponential correlation between bond 
valence and bond length 

vij = exp{{R°jj - Rij)/bij) (2) 

where Rjj is the observed bond length while R^j 
and bij are empirically fitted bond valence parameters 
for I-J bond, bjj measures the bond softness and 
is usually taken as a universal constant of 0.37A, 
while i?5,/ is experimentally determined from structural 
data of related materials[17]. In the present study, 
the disordered alloy offers abundant structural data. 
Therefore R'j j and bij are htted to DFT calculations 
instead. The bond-angle variation depends on the 
higher-order terms of orbital interactions in the 
perturbation expansion. In general the bond bending 
force is weaker than the bond stretching force. In the 
present study, an empirical relation[14] is used for the 
crude prediction of anion-cation-anion angles 

6icj = 109.5 + k{vci + vcj - Vc/2) (3) 

where k is an empirical constant (equal to 15.3° per 
valence unit (v.u.) in Ref. [14]), vci and vcj are 
the bond valences of the two ligand bonds, and Vc is 
the ionic valence of the central cation. Finally, taking 
into account the constraints of bond lengths and bond 
angles, the tetrahedrally coordinated alloy lattice is 
over-constrained. A cost function can be assigned to 
the constraints in order to perform the pre-relaxation. 

3. Computational Methodology 

Here a brief outline is given of the computational meth¬ 
ods. Details are in Ref. [2]. An Ising-type model 


Hamiltonian for the (GaN)i_x(ZnO)x alloy was con¬ 
structed by Li[18] using a DFT-based cluster expan¬ 
sion method[19][20][21]. Monte-Carlo simulations are 
then performed on the constructed cluster expansion 
model using the AT AT package[22][23][24]. For each 
(x, T) of interest, a thermodynamic ensemble of config¬ 
urations is generated. For each configuration the bond 
valences can then be determined using BVM. At this 
stage only the site occupancies are needed. One could 
in principle minimize the measure of the total energy 
E = a J2i j with respect to the set of bond va¬ 
lences {vij}. Unlike the Ising-type cluster expansion 
model, BVM model is essentially short-ranged since 
the set of bond valences {vjj} forms a interactive net¬ 
work. An iterative scheme [25] is used to apply the 
equal valence rule, which generally yields better com¬ 
putational efficiency. Finally the bond-length distribu¬ 
tion and bond-angle variation are obtained using the 
empirical correlations introduced earlier in this section. 

For most of the results presented in this paper, 
the Perdew-Burke-Ernzerhof (PBE)[26] version of the 
exchange-correlation functional is used. Kohn-Sham 
wavefunctions are expanded in a variationally opti¬ 
mized double-^ polarized (DZP) basis set, as imple¬ 
mented in the SIESTA package[27]. Ga-3d and Zn- 
3d electrons are treated explicitly as valence electrons. 
The fc-point mesh is chosen to be equivalent to a 
6x6x4 mesh for the 4-atom wurtzite unit cell. Pseu¬ 
dopotentials for all the atomic species are available 
from the SIESTA homepage[28]. DFT calculations 
are performed for two reasons: (1) The correlations 
Pjj oc vjj, Tlij oc Sij and vjj ^ Sjj are crucial for 
the interpretation of BVM and are therefore exam¬ 
ined first; (2) The BVM parameters are to be fitted to 
DFT calculations, after which bond-length distribution 
and bond-angle variation can be predicted by BVM. I 
construct three representative 432-atom supercells at 
X = 0.25, 0.5 and 0.75 for the former purpose, and 
use a thermodynamic ensemble equilibrated at the ex¬ 
perimental synthesis temperature T = 1,123K[1] with 
72-atom supercells for the latter purpose. 

4. Results and Discussions 

4-T Examination of BVM 

In Fig. 1 the correlations Pjj cx vjj for different types 
of bonds are shown. One should keep in mind that the 
Mulliken population Pjj has no strict physical sense 
due to its sensitivity to the atomic basis set used in 
the projection. Therefore in present study, only the 
qualitative correlation is discussed. The correlation 
Hij oc Sij is in reality adopted in the extended Hiickel 
method [29] where the off-diagonal Hamiltonian matrix 
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Table 1: Bond valence parameters. 













GaN 

original BVM[17] 
GaO ZnN 

ZnO 

GaN 

fitted to DFT 
GaO ZnN 

ZnO 

pO 

bij 

:(^) 

(^) 

1.84 

1.73 1.77 

0.37 

1.704 

1.844 

0.357 

1.755 

0.391 

1.831 

0.268 

1.756 

0.312 


(A) 

1.946 

- 

1.960 

1.947 

- 

- 

1.972 

Expt. 

(A) [32] 

1.95 

- 

1.977 

1.95 

- 

- 

1.977 











elements Hij are approximated by the corresponding 
diagonal Hamiltonian matrix elements and the overlap 
integral through Hij = KSijiHu + Hjj)/2. In Fig. 2 
the correlations Hij oc Sij between the first C numerical 
atomic orbitals of different species are shown. Since 
Ga and O are more electronegative than Zn and 
N respectively, Ga-4s and 0-2p have deeper atomic 
energy levels than Zn-4s and N-2p. This explains why 
the Ga-0 curve lies higher than the Zn-N curve. In 
Fig. 3 the correlations vij ~ Sij for different types 
of bonds are shown. The linearity of the correlations 
validates the interpretation of BVM proposed in the 
present study. 



0.00 0.25 0.50 0.75 1.00 1.25 

V„ 


Figure 1: The correlation P/j oc vij for different types 
of bonds. 

The ability of BVM to predict the bond-length 
distribution relies significantly on the empirical 
correlation vjj = exp {{R'jj — Rij) /bij), the quality 
of which should be examined first. To yield accurate 
structural properties, DFT calculations are performed 
using the QUANTUM ESPRESSO package [30] with 
the PBEsol functional [31]. The lattice constants of 
GaN and ZnO are well reproduced [2]. In Table 1, 
the original (tabulated in Ref. [17]) and fitted-to-DFT 
bond valence parameters are listed. As a sanity check, 
bond lengths of compound GaN and ZnO (labeled 
as calculated with the two sets of bond valence 



Figure 2: The correlation Hij oc Sij between the first 
numerical atomic orbitals of different species. 



Figure 3: The correlation vij ~ Sjj for different types 
of bonds. Sij is dehned as 


parameters are also listed. As the fitting procedure 
releases the freedom of the bond softness bjj, an 
overall improvement is observed for the fitted-to-DFT 
set of bond valence parameters. Fig. 4 shows the 
correlation between the DFT-calculated bond lengths 
and the BVM-predicted bond valences. Bond-length 
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distribution is predicted by BVM with good accuracy. 
The prediction of bond-angle variation is less accurate, 
as shown in Fig. 5. The fitted bond valence parameters 
k for Ga and Zn are 18.1°/v.u. and 20.1°/v.u. 
respectively. In order to perform pre-relaxation, one 
can simply add a penalty function to bond-length 
distribution and bond-angle variation. A large penalty 
to bond-length distribution is suggested while bond 
angles are subject to change. 



Bond valence 


Figure 4: Correlations between the DFT-calculated 
bond lengths and the BVM-predicted bond valences. 
The solid red lines represent the fitted correlations. In 
each figure the number of data points drawn is reduced 
by a factor of ten. 



Figure 5: Correlations between the DFT-calculated 
bond angles and the BVM-predicted bond valences. 
The solid red lines represent the fitted correlations. In 
each figure the number of data points drawn is reduced 
by a factor of ten. 


Jj..2. Bond-length distribution and bond-angle 
variation 

Ref. [2] shows that there is strong SRO in 
the (GaN)i_x(ZnO)x alloy. Although a completely 
random alloy may not be achievable under common 
experimental conditions, the degree of randomness 
introduced is influenced by the experimental methods 
adopted in growing the samples. For example. 


the kinetics of mixing may be facilitated by high- 
pressure [33]. To represent the thermodynamic 

ensembles at {x,T), Monte-Carlo simulations are 
performed on a DFT-based cluster expansion model. 
Given the site occupancies of each configuration, 
bond-length distribution and bond-angle variation are 
then predicted using the fitted-to-DFT bond valence 
parameters. The temperature measures the degree 
of randomness. Temperatures of 1123K, 2000K, 
5000K and 20000K represent short-range ordered 
(SRO), disordered (DISl and DIS2) and random 
(RAN) alloy respectively. Fig. 6 shows the bond- 
length distributions at various temperatures. As 
the temperature is raised, the peak of bond-length 
distribution shifts slightly in the direction of shorter 
bond length. The shift of the peak position is small, 
and can be easily overwhelmed by other factors such 
as thermal expansion, which is not considered here. In 
the meanwhile the width of bond-length distribution 
becomes broader with increasing randomness. In 
Fig. 7, bond-length distributions of different types 
of bonds are shown. Upon mixing, the Ga-N bond 
shrinks while the Zn-0 bond expands relative to the 
bond lengths in the corresponding compounds. From 
SRO alloy to RAN alloy, the shift is toward shorter 
bond length for Ga-N, barely temperature-dependent 
for Ga-0 and Zn-N, and is reversed to the longer 
bond-length direction for Zn-0. This unusual tendency 
of bond-length distribution is a consequence of the 
non-isovalent nature of the alloy, and can be easily 
interpreted in terms of bond valence. One consequence 
of elevating the degree of randomness is to enhance 
the statistical presence of the energetically unfavored 
valence-mismatched Ga-0 and Zn-N pairs. In the 
language of BVM, for a cation-anion pair, enhancing 
the presence of N(0) neighbors around the cation and 
Ga(Zn) neighbors around the anion will drain(pour) 
bond valence from(into) the cation-anion pair and as a 
result the bond is lengthened(shortened). Of particular 
importance is the Zn-N bond-length distribution due 
to the decisive role of Zn3d-N2p repulsion on the top 
of the valence band. In Ref. [2], an almost linear 
band gap reduction upon increasing the ZnO content 
for the short-range ordered alloy is observed. Since 
the p-d repulsion is inversely proportional to the bond 
length, upon increasing the ZnO content a shortened 
Zn-N bond-length distribution is expected, which is 
confirmed by the BVM prediction shown in Fig. 8. 
The stronger p-d repulsion pushes the top of the valence 
band, resulting in the linear band gap reduction. Fig. 
9 shows the anion-cation-anion bond-angle variation 
of short-range ordered alloy at x = 0.5. The N-Ga-N 
angle expands while the O-Zn-0 angle shrinks relative 
to the ideal tetrahedral angle 109.5°, which can be 
explained by noticing that the bond valence of the 
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ligand cation-0 bond is generally smaller than that of 
the ligand cation-N bond. For Fig. 7-9, see Ref. [2] for 
the DFT-calculated more reliable but less statistical 
predictions. 



■ ^ I ■ I , i^ **^*^ - *.i; 

1.7 1.8 1.9 2.0 2.1 2.2 2.3 


Bond length (A) 

Figure 6: Temperature dependence of bond-length 
distribution at a; = 0.5. 



1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 

Bond length (A) 


Figure 7: Bond-length distributions of short-range 
ordered alloy and random alloy at a; = 0.5. The 
vertical dotted lines mark the bond lengths of the 
corresponding compounds. 


Energetics 

As for the energetics, DFT total energy calculations 
are performed on 170 structures selected from the 
r = 1,123K thermodynamic ensemble over the full 
range of composition. The formation energies are also 
calculated using the valence loop rule (i.e. the mea¬ 
sure of energy E = a ..j 'I’h)- The results are shown 
in Fig. 10. The fitted parameter a = 1.07 is con¬ 
sistent with that of Ref. [9]. The power of BVM is 



Bond length (A) 

Figure 8: Zn-N bond-length distribution at various 
ZnO content. 



Figure 9: The anion-cation-anion bond-angle variation 
of short-range ordered alloy at a; = 0.5. 

shown by the accurate reproduction of the energetics. 
However, BVM fails to reproduce the ordered superlat¬ 
tice (GaN)i(ZnO)i ground state at a; = 0.5[18], possi¬ 
bly due to the nearest-neighbor short-range nature of 
BVM itself, i.e. the wurtzite and zincblende structures 
are indistinguishable from one another in BVM. The 
formation energy of (GaN)i(ZnO)i predicted by BVM 
is positive, while that predicted by DFT is slightly 
negative[18]. The discrepancy should not affect any 
conclusion drawn in present study since only the dis¬ 
ordered phase is concerned. 

Inclusion of vibrational entropy in the first-principles 
alloy phase diagram calculation is a long-standing chal¬ 
lenge. The main difficulty lies in the conflict between 
the requirement for a large supercell and the expen¬ 
sive computational cost associated with it. The prob¬ 
lem is partly alleviated by the SQS approach[34][35][2]. 
Another idea is to use a bond-length-dependent trans¬ 
ferable force constant [11] [36], where the bond stiff¬ 
ness is predicted from the bond length and the chem¬ 
ical identity of the bond. The present study reveals 
the strong correlation between bond valence and bond 
length, which indicates the possibility of using bond 
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Figure 10: Comparison of DFT-calculated formation 
energies with BVM-predicted formation energies. 



Bond length (A) 


< 


Figure 11: Correlation between stretching bond 
stiffness and bond length. is shown in black, and 
is shown in red. 


valence instead of bond length as a predictor for bond 
stiffness. Such extension will release the estimation 
of nearest-neighbor force constants from the require¬ 
ment for the knowledge of the relaxed geometry of 
a conhguration. Fig. 11 shows the dependence of 
stretching bond stiffness and cpp^^ on bond length, 
where a refers to the bond-stretching direction and 
/, J are nearest neighbors. The bond stiffness cal¬ 
culations are performed on selected 72-atom super¬ 
cells with a displacement of 0.02A from the relaxed 
atomic coordinates along each bond direction. While 
a linear bond stiffness vs bond length relationship is 
suggested in bond-length-dependent transferable force 
constant approach [11] [36], an exponential correlation 
(similar with that between bond valence and bond 
length) seems to fit better according to the present 
study, which is consistent with the interpretation that 
bond valence measures bond strength. Bond stiffness 
depends on bond length in a similar manner. The most 
covalent Ga-N bond is the stiffest, while the most ionic 
Zn-0 bond is the softest. 

For isovalent III-V semiconductor alloys, the widely 
used Keating valence force field (KVFF) model[37] 
yields generally good accuracy[38][39][40][41][42][43]. 
In KVFF, the force constants are related to the macro¬ 
scopic elastic constants, and therefore can be accu¬ 
rately determined experimentally. Also the isovalent 
nature of IIl-V semiconductor alloys guarantees good 
transferability from compound semiconductors to the 
corresponding alloy. For non-isovalent semiconductor 
alloys, the transferability no longer holds, for the ap¬ 
parent reason that there exists no wurtzite GaO or 
ZnN. The present study offers an alternative way of 
accurately reproducing the energetics of non-isovalent 
semiconductor alloys with BVM, where only site oc¬ 


cupancies are needed. The non-isovalent nature is well 
captured by the valence sum rule. An extension of the 
BVM leads to the modelling of an atomistic potential. 
In present study, the relaxation energy is assumed to 
consist of three parts: 

^relax —O^ ^ ^ ^IJ T ^ ^ kl ^ ^ ^o) 

I.J I=Ga,Zn J 1 .J 2 

+ X! ~ ^OJ 

I=N.O V J 

The first term is simply the valence loop rule, and 
the second term is the harmonic angle potential. 
The third term accounts for large lattice relaxations 
by penalizing deviations from the bond valence 
conservation and is important for reliable molecular 
dynamics simulations [5] [6]. In the fitting procedure 
each relaxed structure is expanded and contracted by 
1%. Fitting parameters A:n,o and /Jca.Zn are found to 
be negligible. In Figure 12 the comparison between 
DFT-calculated and BVM-fitted formation energies 
is shown. The agreement is generally satisfactory. 
Further studies will involve refinement of the atomic 
potential. 

5. Discussion and Conclusions 

A physical interpretation of BVM is discussed from the 
computational perspective. The underlying assump¬ 
tions and correlations within BVM are revealed by 
DFT calculations on the non-isovalent semiconductor 
alloy (GaN)i_x(ZnO)x. Bond-length distribution and 
bond-angle variation are predicted by fitting BVM em¬ 
pirical relations to reliable DFT-calculated structural 
data. The unusual relaxations associated with the non- 
isovalent nature of the alloy are explained. Effects of 
SRO on bond-length distribution and bond-angle vari¬ 
ation are also discussed. The energetics is accurately 
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Figure 12: Comparison between DFT-calculated and 
BVM-fitted formation energies. 

reproduced by BVM. The connection between bond va¬ 
lence and stretching bond-length-dependent transfer¬ 
able force constant is revealed. A tentative improved 
bond valence potential is proposed. In principle, the 
methods of the present study should also be applicable 
for other non-isovalent semiconductor alloys. 
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